#整體設定，含載入套件
source("https://raw.githubusercontent.com/ChungPingCheng/R4BS2/main/R4BS_setup.R")

#資料來自於 TIMSS 2019 年台灣資料，納入變項無缺漏值者
#讀檔案
dta <- read.table(file="../Data/TIMSS2019TW.txt", 
                  header = TRUE, stringsAsFactors = TRUE)

#看資料結構與前六筆（類別變項與連續變項分開呈現）
#程式報表4.1
str(dta)

#排除無關變項（有「科學「兩字的）
#把父母教育各個水準順序定下來
dta <- dta |>   
  dplyr::select_at(vars(-contains("科學"))) |>
  dplyr::mutate(父母教育 = factor(父母教育, 
    levels = c('國小以下', '初中', '高中','專科','大學以上')))

#呈現類別變項前六筆
#程式報表4.2前
dta |>  dplyr::select_if(is.factor) |> head()

#呈現連續變項前六筆
#程式報表4.2後
dta |>  dplyr::select_if(is.numeric) |> head()

#看資料基本統計
#程式報表4.3
dta |> 
  gtsummary::tbl_summary(by = '父母教育',
       statistic = list(all_continuous() ~ "{mean} ({sd})")) 

#圖示不同父母教育程度下的數學分數平均數，加上信賴區間
#圖4.1
ggplot(data = dta, 
       aes(x = 父母教育, 
           y = 數學)) +
  stat_summary(fun.data = 'mean_cl_boot') +
  scale_y_continuous(breaks = seq(500, 660, by = 20)) +
  geom_hline(yintercept = mean(dta$數學) , 
             linetype = 'dotted') +
  coord_flip() +
  labs(x = '父母教育程度', y = '數學平均分數',
       title='數學分數與父母教育程度') 

#父母教育的效果或許會是教育資源造成的，畫圖看看
#圖4.2
ggplot(data = dta, aes(y = 數學, x = 教育資源)) +
  geom_point(pch=1, size=rel(.3), alpha=.3)+
  stat_smooth(method = 'lm', 
              se = F, 
              formula = y ~ x,
              col=1,
              linewidth=.8) +
  facet_wrap(vars(父母教育), nrow=1)+
  labs(x = '教育資源', y = '數學分數')

#以 tidy 方式進行分析並輸出報表
#程式報表4.4
dta |> 
  lm(formula = 數學 ~ 父母教育, data = _) |> 
  broom::tidy(conf.int = TRUE) |>
  insight::format_table(ci_brackets = c("(", ")")) |>
  export_table(format = "markdown", caption = "", align = "lcccrr")

#這是原始輸出
#程式報表4.5
dta |> 
  lm(formula = 數學 ~ 父母教育, data = _) |> 
  summary()

#ANOVA 式報表
#程式報表4.6前
anova(m1 <- lm(數學 ~ 父母教育, data = dta))

#把教育資源加進模型
#程式報表4.6中
anova(m2 <- update(m1, . ~ . + 教育資源, data = dta))

#或許不是父母教育而是教育資源造成，這邊只考慮教育資源
#程式報表4.6後
anova(m3 <- update(m2, . ~ . - 父母教育,  data = dta))

#比較m3,m2
anova(m3, m2)

#比較m1,m2
anova(m1, m2)

#將結果放在一個list中，等一下比較方便抓結果
mod_r2 <- sapply(list(m1, m2, m3), broom::glance) |>
  as.data.frame() |>
  dplyr::slice(1) |> 
  t() |> 
  unlist()
mod_r2

#比較在控制教育資源下，父母教育的效果
#程式報表4.7 前半
(mod_r2[2] - mod_r2[3])/mod_r2[2]

#比較在控制父母教育下，教育資源的效果
#程式報表4.7 後半
(mod_r2[2] - mod_r2[1])/mod_r2[2]

#輸出最終模型結果
#程式報表4.8
dta |> 
  lm(formula = 數學 ~ 父母教育 + 教育資源, data = _) |> 
  jtools::summ(confint = TRUE)

#以 ggplot 畫效果
#將截距去除，畫更易懂起來
#圖4.3
tidy(m2, conf.int = TRUE) |> 
  dplyr::slice(-1) |>
ggplot() + 
  aes(estimate, term, xmin = conf.low, xmax = conf.high, height = 0) +
  geom_point(size=rel(3)) +
  geom_vline(xintercept = 0, linetype="dashed") +
  geom_errorbarh() +
  labs(x = '估計值和信賴區間', 
       y = '迴歸模型變項(去除截距)', 
       subtitle = '依變項:數學分數')

#以 plot_summs 畫效果
#將截距去除，畫更易懂起來
#圖4.4
jtools::plot_summs(m2) +
  labs(x = '估計值和信賴區間', 
       y = '迴歸模型變項(去除截距)', 
       subtitle = '依變項:數學分數')

#預備要畫預測時，先看一下教育資源數值差異會到多少
#程式報表4.9
with(dta, 教育資源) |> quantile()

#畫預測值，移除 colors='bw' 可以看到彩色圖
#圖4.5
sjPlot::plot_model(m2, type='pred', colors='bw',
  terms=c("父母教育","教育資源[9.036, 11.487]")) +
  labs(y="數學分數估測值",
       title="數學 ~ 父母教育 + 教育資源")

#利用 gglm 套件作模型診斷 
#圖4.6
gglm::gglm(m2, theme = ggplot2::theme_minimal())

#利用 performance 套件作模型診斷 
#圖4.7
performance::check_model(m2)

#殘差圖
#畫圖前用 broom 把迴歸分析的預測值、殘差與槓桿值放進資料
#圖4.8, 上方左圖
p1 <- broom::augment(m2) |>
ggplot() + 
  aes(.fitted, .std.resid) +
  geom_hline(yintercept = 0, col='gray') +
  geom_point(pch = 1, size=rel(.5), alpha=.5) +
  geom_smooth(method='rlm', 
              formula = y ~ x, 
              se = FALSE,
              linewidth=.5, 
              col='black')+
  facet_wrap(vars(父母教育), ncol=1) +
  labs(x="數學分數估測值", 
       y="標準化殘差",
       title = "數學 ~ 父母教育 + 教育資源")

#看看殘差的 Q-Q 圖，依父母教育。檢視常態假設
#圖4.8, 下方
p2 <- broom::augment(m2) |>
ggplot() +
  aes(sample=.std.resid)+
  facet_wrap(vars(父母教育), nrow=1)+
  stat_qq(pch=1, size=rel(.5), alpha=.5)+
  stat_qq_line(linewidth=.5)+
  labs(x = '常態位數', y = '標準化殘差')

#呈現槓桿值與標準化殘差
#區分父母教育
#圖4.8, 上方右圖
p3 <- broom::augment(m2) |>
ggplot() +
 aes(x = .hat, y = .std.resid) +
 geom_point(pch=1, alpha=.5, size=rel(.5))+
 geom_hline(yintercept = 0, linetype = 'dotted') +
 facet_wrap(vars(父母教育), ncol=1) +
 labs(x = '影響值', y = '標準化殘差')

#三張圖併成一張呈現
(p1 + p3) / p2

#看看影響值
#程式報表4.10
broom::augment(m2) |> dplyr::select(".hat") |> summary() 

#底下要呈現多個連續解釋變項時的情形
#看看個人變項的可能效果，把跟數學有關的部分取出來
#看看基本統計量
dta |> 
 dplyr::select_at(vars(starts_with("數學"))) |>
 dplyr::select_if(is.numeric) |> 
 colMeans()

#呈現兩兩散佈圖
#以corrplot 套件，利用圖形顯示相關大小
#圖4.9
dta |> 
 dplyr::select_at(vars(starts_with("數學"))) |>
 dplyr::select_if(is.numeric) |> cor() %T>%
corrplot(corr = ., method = 'ellipse', order = 'hclust', addrect = 4,
         type = 'upper', tl.pos = 'd') %>%
corrplot(corr = ., add = TRUE, type = 'lower', method = 'number',
         order = 'hclust', col = 'black', 
         diag = FALSE, tl.pos = 'n', cl.pos = 'n')

#放進三個解釋變項
#程式報表4.11前
dta_math <- dta |> 
 dplyr::select_at(vars(starts_with("數學"))) |>
 dplyr::select_if(is.numeric) 
m4 <- lm(數學 ~ 數學興趣 + 數學評價 + 數學信心, data = dta_math)
m4 |> broom::glance() |> dplyr::select(1:6)

#利用 lm.beta套件，計算標準化迴歸係數
#程式報表4.11後
m4 |> lm.beta::lm.beta() |> summary()

#看效果
#圖4.10
jtools::plot_summs(m4) +
  labs(x = '估計值和信賴區間', 
       y = '迴歸模型變項(去除截距)', 
       subtitle = '依變項:數學分數') 

#看效果
#圖4.11
p1 <- jtools::effect_plot(m4, pred = '數學興趣', 
                    interval = TRUE, 
                    point.size = .1, point.alpha = .3,
                    plot.points = TRUE)
p2 <- jtools::effect_plot(m4, pred = '數學評價', 
                    interval = TRUE,  
                    point.size = .1, point.alpha = .3,
                    plot.points = TRUE)
p3 <- jtools::effect_plot(m4, pred = '數學信心', 
                    interval = TRUE, 
                    point.size = .1, point.alpha = .3,
                    plot.points = TRUE)

p1 + p2 + p3

#看看控制數學信心與數學評價後，數學興趣的效果
#程式報表4.12前
summary(m5 <- update(m4, . ~ . - 數學興趣 , data = dta_math))

#程式報表4.12後
anova(m5, m4)

#看看不同模型的預測以及信賴區間
#圖4.12
jtools::plot_summs(m4, m5, model.names = c(formula(m4), formula(m5)))  +
  labs(x = '估計值和信賴區間', 
       y = '迴歸模型變項(去除截距)', 
       subtitle = '依變項:數學分數') +
  theme(legend.position="bottom")